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I.    INTRODUCTION 

A .  Background 

The  Northrop  Corporation,  under  Air  Force  funding,  has  developed 
a  finite  element  digital  computer  code,  called  BR-1,  for  predicting  the 
inelastic,  large  deflection,  transient  response  of  combat  aircraft  skin- 
rib-stringer  structures  when  subjected  to  internal  air  blast  loading. 
The  finite  elements  considered  are  flat  rectangular  plates  and  beam 
stiffeners.  The  theory,  user's  manual  and  code  listing  are  given  in 
References  1  and  2.  The  Air  Force  Flight  Dynamics  Laboratory  wanted 
the  BR-1  code  modified  so  that  it  could  be  used  to  predict  the  response 
of  aircraft  fuel  tank  walls  when  subjected  to  fluid  pressures  due  to 
projectiles  passing  through  the  fuel  in  the  tank.  The  intense  pressure 
and  momentum  in  the  fuel  due  to  the  penetrating  projectile  is  referred 
to  as  the  hydraulic  ram  loading.  This  report  describes  the  modifica- 
tions to  the  IBM  version  of  the  BR-1  code  to  account  for  the  fluid 
(fuel)  -  structure  (tank  wall)  interaction  that  occurs  when  bullets  and 
metal  fragments  penetrate  into  aircraft  fuel  tanks.  The  modified  code 
is  called  BR-1HR.  The  interaction  between  the  compressible  fluid  and 
the  structure  is  approximated  by  the  piston  theory.  The  code  can  also 
be  used  for  many  other  compressible  fluid-structure  interaction  problems 

B.  Piston  Theory 

The  total  nonlinear  problem  of  the  response  of  a  tank  containing  a 
fluid  and  subjected  to  a  high  speed  penetrating  projectile  is  extremely 
complex  and  presently  defies  analytical  treatment.   In  general,  the 


equations  for  the  fluid  stresses  and  motion  are  coupled  to  those  for  the 
wall  stresses  and  motion  due  to  the  continuity  at  the  fluid-structure 
interface  (3).  One  procedure  for  approximating  the  fluid-structure  in- 
teraction phenomenon  is  the  piston  theory  (h) .  This  theory  has  been  in 
use  since  the  early  19^0' s  when  it  was  applied  to  the  study  of  the 
effect  of  underwater  explosions  on  ship  plates.   It  provides  the  correct 
solution  to  the  one -dimensional  propagation  of  stresses  in  an  acoustic 
medium  due  to  a  moving  boundary.  Several  recent  studies  have  been  made 
to  determine  its  accuracy  when  applied  to  two  dimensional  fluid-structure 
interaction  problems  (*+,5). 

Application  of  the  piston  theory  to  the  interaction  problem  allows 
the  structure  equations  and  fluid  equations  to  be  uncoupled.  The  response 
of  the  wall  is  computed  using  the  conventional  structural  response  equa- 
tions, with  the  normal  pressure  on  the  wall  p  given  by 

p  =  Po  +  pc  (v±   -  w)  (1) 

where  p  and  v  are  the  incident  pressure  and  velocity  of  the  fluid  at 
o      1 

the  wall  respectively,  P  is  the  fluid  density,  c  is  the  acoustic  velocity 
in  the  fluid,  and  w  is  the  wall  velocity*.  The  pressure,  p  ,  and 
velocity,  v. ,  are  the  pressure  and  velocity  that  would  exist  in  the 

fluid  if  the  interface  was  not  there,  i.e.,  p  and  v.  do  not  contain 

'        o      1 

any  "local"  reflected  effects.  However,  effects  on  p  and  v.  due  to 

o      1 

earlier  reflections  from  other  walls  and  free  surfaces  should  be  con- 
sidered. In  other  words,  p  and  v.  are  the  loading  components  due  to 


*A  dot  above  a  variable  denotes  a  derivative  with  respect  to  time 


the  free  field  and  the  scattered  effects.  The  loading  component  due  to 

the  wall  velocity  w  is  called  the  radiation  pressure. 

C.  The  NWC  Hydraulic  Ram  Computer  Code 

In  order  to  use  the  piston  theory  to  compute  the  tank  wall  response, 

it  is  necessary  to  know  the  incident  fluid  pressure  p  and  velocity  v 

o  i 

over  the  entire  fluid-wall  interface  as  a  function  of  time.   In  conjunc- 
tion with  this  project  Lundstrom,  at  the  Naval  Weapons  Center,  has 
developed  a  digital  computer  code  that  predicts  the  fluid  pressures  and 
velocities  p  and  v.  throughout  a  rectangular  body  of  fluid  due  to  a 
penetrating  ballistic  projectile.  The  model  is  based  upon  replacing 
the  projectile  by  a  line  of  sources  whose  strength  is  determined  by  an 
energy  balance  between  the  kinetic  and  potential  energy  of  the  fluid  and 
the  energy  loss  due  to  drag  forces  on  the  projectile.  Reflections  from 
the  structure -fluid  interface  are  accounted  for  by  considering  the  fluid 
boundary  to  be  either  stress  free  or  rigid*.   An  extensive  series  of 
tests  were  performed  at  the  Naval  Weapons  Center  to  obtain  detailed 
pressure  measurements  for  a  variety  of  projectiles  under  a  wide  range  of 
impact  conditions.  This  data  allowed  the  selection  of  important  para- 
meters such  as  tumbling  distance,  jacket  stripping,  etc.,  to  be  entered 
into  the  code.  A  description  of  the  code  and  the  instructions  for 
operation  are  given  in  Reference  7.  This  code  provides  the  values  for 

p  and  v.  at  user  specified  locations  over  the  structure -fluid  interface 

o      1         * 

for  the  time  span  of  interest. 


*-  A  study  of  the  one-dimensional  reflection  of  step  pressure  waves  from 
typical  aircraft  fuel  tank  walls  indicates  that  the  stress  free  surface 
provides  the  more  accurate  approximation  (6) 


II.  MODIFICATION  OF  THE  BR-1  CODE 

A.   Incorporation  of  the  Piston  Theory 

The  BR-1  code  has  an  option  for  the  user  to  input  a  time  varying 
pressure  on  each  panel  element.   In  the  piston  theory  this  pressure  is 
the  p  +  pcv.  of  Eg.  1.  The  other  contributor  to  the  wall  loading  given 
"by  Eq.  1  is  w,  the  wall  velocity.  Since  the  BR-1  code  does  not  include 
damping  effects,  it  is  necessary  to  add  the  damping  term  pew  to  the 
equations  of  motion. 

The  BR-1  code  solves  the  set  of  equations  (Eqs.  1  and  2,  Ref.  l) 

[M]  {q*}  =  {F}  -  {P}  -  [H]  {q*}  =  {C}  (2) 

for  the  vector  of  global  nodal  generalized  displacements  {q*}  as  a  function 
of  time.  These  generalized  displacements  define  the  motion  of  the  walls. 
The  vector  {F}  consists  of  global  generalized  external  and  body  forces 
at  the  nodes  of  the  elements.  The  matrix  [M]  is  the  mass  matrix. 

The  wall  pressure  p  given  by  Eq.  1  causes  external  forces  at  the 
nodes.  The  external  generalized  forces  at  the  nodes  of  each  element  in 
the  local  coordinate  system,  {f},  is  given  by  (Eq.  A-k7 ,   Ref.  l) 


«-LWW« 


(Aout ) 
T 
where  L  BJ   is  the  transpose  of  the  matrix  of  shape  functions  [N] , 

evaluated  at  the  surface  of  the  element,  Aout  is  the  surface  of  the  ele- 
ment, and  {T  }  is  the  vector  of  applied  surface  tractions  and  moments. 


The  order  of  (T  ]  is  a  5x1  vector.  Due  to  the  fluid  pressure  loading 


N  H 

<  T„  >  =  <  0  > 


I  3. 


vP^ 


00 


where  the  subscripts  denote  the  coordinates,  (T  is  normal  to  the  ele- 
ment), and  p  is  given  by  Eq.  1.  The  fourth  and  fifth  elements  of  [T  } 
correspond  to  applied  moments  per  unit  middle  surface  area  of  the 
element,  and  are  zero  here.  The  numerical  intergration  of  Eq.  3  can 
be  accomplished  by  Gaussian  quadrature.  However,  the  [f]  is  obtained 
in  the  BR-1  code  in  a  more  approximate  way  by  using  a  lumping  approach 
at  the  nodes  of  the  element,  as  is  done  for  the  mass  matrix  evaluation, 
in  order  to  save  computation  time.  Thus,  according  to  Eq.  B-92  of 
Ref .  1,  the  external  force  vector  at  the  rth  node  of  the  nth  element 
is  given  by 

nr  I  .. 

0 

0  /nr 

where  p  is  the  magnitude  of  the  pressure  on  the  element*,  and  D 

vn  r-  '      nr 


V'    2  2 
(l+0n-+9„)    where  fin  and  0^  are  the  fourth  and  fifth  local  general- 
^   ul  w2  nr       1      2 

ized  displacements  at  the  rth  node.  They  appear  here  because  the 
pressure  is  defined  in  BR-1  as  the  pressure  normal  to  the  deformed 
surface.  The  quantities  (x  -x  )  and  (y  -yp)  are  the  dimensions  of  the 
rectangular  element. 

The  pressure  p  in  the  piston  theory  is  given  by  Eq.  1,  i.e. 

P  =  P   +  (pc)  v.   -  (pc)  w  (6) 

rn    on    F  n  in    K  n  nr 


*  The  assumption  is  made  in  the  programming  of  BR-1  that  the  pressure 
is  uniform  over  each  element.  This  is  contrary  to  the  theoretical 
presentation  where  the  pressure  is  defined  at  each  node  point . 


where 

|  pC  if  the  nth  element  is  in  contact  with  the  fluid 

'n  =  J 

(  0  if  the  nth  element  is  not  in  contact  with  the  fluid. 


(PcL  = 


The  variables  p   and  v.   can  be  determined  from  the  NWC  computer  code 
on      in 

for  each  element  for  the  time  span  of  interest  prior  to  the  computation 

of  the  wall  response.  This  data  is  then  input  as  the  known  external 

pressure.  The  variable  w   is  an  unknown  dependent  variable  and  is  part 

of  {q*}.  Hence,  it  must  be  incorporated  into  the  equations  of  motion, 

Eq.  2. 

The  ff  }  due  to  w   is  given  by 
1  nrJ        nr 


■M  H, 


,92 
/pCx   w    1-6 
ff    =   \     /      \  /"  V     )n     nr  I      1 

1  nrJ       ""^T  D       1 

'  0 
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(7a) 


nr 


The  global  force  vector  at  the  rth  node,  [F]  ,  is  related  to  [f  }  in 
the  form 

[F]  =  E   [J  ]  T  ff  1  (7b) 

where  [J  ]  is  the  transformation  matrix   from  the  global  coordinate 
system  to  the  nth  element  local  coordinate  system,  v  means  a  summation 
over  all  elements  containing  the  node  r,  and  oi  < — »  r  means  node  a   corre- 
sponds physically  to  node  r.  Since  the  .wall  velocity  in  Eq.  7a  is  given 
in  terms  of  the  local  coordinates,  it  must  be  converted  to  global 
coordinates.  Thus,  according  to  Eqs .  A-77  and  B-3  of  Reference  1, 

w   =  LJ3J  fq*}  (8) 

nr   L  nJ  L   Jr 


where  I  J  I  is  the  third  row  of  fj  "I.  Thus,  the  global  external  force 
L  nJ  L  nJ         3 

vector  at  the  rth  node  [F]  becomes 


{F}r  =  'I  g   (pOn  [j/  (^-XD/g-yl)°  [J3]  i4.j 
"r  nr 


f  ft  ^ 

"ei 
i    ' 

0 

Lo  J 

(9) 


nr 

•      •  • 

Note  that  [F]   is  nonlinear  since  q  and  ft  are  part  of  {q*}  •   If  the 

rotations  0  and  ft  are  neglected  in  the  computation  of  [F]  ,  i.e.,  if 
the  pressure  is  not  truly  normal  to  the  deformed  surface,  [F]  for  the 
total  R  nodes  of  the  structure  can  be  given  in  the  form 


[F]  =  -  [D]  {4*} 


(10a) 


where 


[D]  = 


[D], 


P>! 


and  [D]   is  a  6x6  matrix  given  by 


(10b) 


•[D] 


R 


'X^    x 


yo  y^ 


[d]     =iz    (pc)     ("2-"l)     (^2-^1)     [J  ]x  LJ  Ji 
Ljr        4^      \p'n\  /  /nLnJ      LnJ 


(10c) 
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B.  Method  of  Solution 

The  BR-1  code  solves  for  {q*}  at  discrete  points  in  time  using  the 
explicit  finite  difference  scheme  (Eq.  A-109,  Ref .  l) 

[Aq*}t    =  {Aq*}t   +  At  [qj}   +  (At)2  (q*}  (ll) 

i+1        i  i  i 

where  At  is  the  time  interval  between  two  time  points,  i.e. 

At  =  *1+1  "  *1 

and  {Aq*}t        =  [q*J  "   {l*}t  ^12) 

i+1       i+1        i 

The  acceleration  fq*}   is  obtained  from  Eq.  (2)  in  the  form 

i 

[q*}t   ==  CM]"1  [C}t  (13) 

i  i 

The  {q£}  are  due  to  impulsive  loads  which  are  known  in  the  blast  loading 

problem.   In  the  BR-1  code  {F},  {P},  {q*}  and  {q*}  are  known  at  time  t.. 

Hence,  {Aq*}     and  {q*},     can  be  determined  using  Eqs.  11-13 .  For 
ti+l        ti+l 

our  situation,  {F},   contains  {q*}.  ,   which  is  unknown.  We  could 

i  i 

approximate  {q*}.   with  the  backward  finite  difference  form 


1 


{q*}   =  {Aq*}   /At  (ik) 


1 


If  we  do,  then  {q*}.   becomes  known  at  t.,  {F},  and  hence  [C],  can  be 

i 

determined  at  t . ,  and  the  procedure  used  in  BR-1  is  directly  applicable . 

On  the  other  hand,  if  we  express  {q*}.   in  the  central  finite  difference 

i 
form 


'*\  -  ( 


[q*}t   =    Aq*}t>    -  {Aq*}.      /  (2At)  (l5) 


then  {F}   ,  and  hence  {C}   ,  depends  upon  {Aq}    .   Consequently 

i  i  1+1 

{Aq*}     appears  on  both  the  left  and  right  hand  side  of  Eq.  11.  This 

i+l 
requires  a  new  solution  procedure.  A  detailed  study  of  the  accuracy 

and  numerical  stability  of  these  two  approaches,  and  a  third  approach, 
when  applied  to  a  single  degree  of  freedom,  damped  oscillator  is  pre- 
sented in  the  Appendix.  The  approach  where  {q*}  is  given  by  the  central 
difference  expression,  Eq.  15,  is  the  one  selected  based  upon  the 
accuracy  and  stability  properties  of  this  scheme.   Its  shown  in  the 
Appendix  that  the  maximum  value  of  At  for  a  stable  solution  is  2/^,  where 
^  is  the  highest  natural  undamped  frequency.  This  is  identical  to  the 
stability  limit  on  the  BR-1  procedure. 

Introducting  that  part  of  {F}  due  to  w  given  by  Eq.  10a  into  Eq.  2 
results  in  the  modified  equations  of  motion 

[M]  {q*}  +  [D]  {4*}  =  {C}  (16) 

Replacing  {q*}  with  the  conventional  central  difference  approximation 

i   \     i+l        i       i+l// 
is  equivalent  to  obtaining  [4*}+  from  Eq.  11  with  {q*-}  not  considered,  i.e, 

{q*}t_  =  ({Aq*}t>+   "  {Aq*}t>  j/(At)2  (17b) 

according  to  Eq.  12.  Substituting  Eq.  15  for  {q*}+  and  Eq.  17b  for 

i 

{el*}.!,      into  Eq.    16  leads  to 
i 

{Aq*}t  =   [M  +  D    (At/2)]"1    (    [M  -  D    (At/2)]    {Ad*}        +   {C}t     )(l8) 

i+l  \  i  1' 


10 


which  is  equivalent  to  {Aq*}     given  by  Eqs .  11  and  13  when  D  is  a 


i+1 


null  matrix  and  {q*}  is  not  considered. 

The  mass  matrix  [M]  is  developed  in  BR-1  using  the  lumped  mass 
approach  and  is  given  by  (Eq.  A-97>  Ref .  l) 


m 


[Mj.1 


[Mo] 


0 


0 


py 


(19a) 


where   [M]      is  a  6x6  matrix  given  by 


[Ml      =  S     r J  ]        [m     ]  TJ  1 
L   Jr       n         n  ruJ 


la 


(19b) 


and  [m    is  a  diagonal  matrix  of  the  lumped  mass  at  the  a   node  of  the 

nth  element.   Comparing  Eq.  10b  with  Eq.  19a  reveals  that  the  two 
matrices  [M  +  D  (At/2)]  and  [M  -  D  (At/2)]  occupy  the  same  nonzero  6x6 
locations  as  the  original  matrix  M.  Thus,  the  same  procedure  used  in 
BR-1  to  compute  [M]    can  be  used  to  compute  [M  +  D  (At/2)]"  .   Its 
only  necessary  to  modify  the  elements  of  [M]  by  the  addition  of  the 

damping  matrix  [D]   given  by  Eq.  10c.  The  other  necessary  change  is 

the  addition  of  the  Matrix  [M  -  D  (At/2)]  as  a  product  with  {Aq-*}   . 

i 

Thus  since 
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IT1- 


r  ^ 

-i 

[M2] 

c 

0      ''. 

_ 

''[MR]"1  J 

(20) 


Eq.  18  can  be  expressed  in  the  form 


{Aq*}rt    =  [M  +  Dr  (At/2)]"1  ([1^  -  ^    (At/2)]  fAq*}rt 
i+1  \  i 


(21) 


+  {C} 


rt./  ■ 


r  =  1,  2,...R 
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C.   Program  Changes  and  Modification  Logic 

The  following  routines  of  the  IBM  version  of  BR-1  have  "been  modified 
for  BR-1HR:   MAIN,  MEMBER,  MTERM,  QPLATE,  ST0RE,  DELTT ,  DEFLX,  REVIV1 
and  REVIV2.  Two  new  subroutines  were  created:  DPMASS  and  ADDAMP.   The 
core  size  was  increased  from  250k  bytes  to  290k  bytes. 

The  flow  of  the  logic  of  the  modifications  is  as  follows: 

1.  Compute  [M]   in         ST0RE         (no  change) 

2.  Compute  [D]   in         ST0RE 

3.  Compute  [M]   in         MTERM         (no  change) 

k.      Compute  maximum  time  interval  for  numerical  stability  in  DELTT 
based  on  [M]  (no  change) 

5.  Take  the  inverse  of  [M]  "  to  get  [M]   in  DPMASS  using  UWS 

6.  Compute  [M  +  D   (At/2)]  and  [M  -  D   (At/2)]  in  DPMASS 

7.  Compute  [M  +  D   (At/2)]"1  in  DPMASS  using  INVS 

8.  Compute  [M  +  D   (At/2)]"1  [M  -  D   (At/2)]  in  DPMASS 

9.  Compute  [M  +  D   (At/A)]  "1  [M  -  D   (At/A)]  {Aq*}    in  ADDAMP 

r    r  r    r  "C . 

i 

10.   Compute  rAq*},    using  Eq.  21  in  DEFLX  (no  change) 

in 

The  phrase  "no  change"  means  that  the  original  procedure  was  used. 

When  no  damping  is  considered  the  modifications  and  additions  are 

bypassed. 
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III.   USER'S  INSTRUCTIONS  FOR  BR-1HR 

The  instructions  for  the  use  of  the  BR-1  code  are  given  in  Ref.  2. 
All  of  the  instructions  contained  in  that  volume  also  apply  to  the  modi- 
fied program  BR-1HR.  The  time  step  for  numerical  stability  of  BR-1HR 
is  identical  to  that  of  BR-1.  The  additional  instructions  required  to 
use  BR-Uffi  are  as  follows: 

1.  Problem  Control  Card  (page  k,   Ref.  2) 

IHR  (1 5,  Col.  66-70)  -  IHR  =  0,  no  fluid  is  involved;  the 

original  BR-1  code  is  used.   IHR  =  1, 
(follows  IREV) 

fluid  is  involved,  the  modifications 

are  used. 

2.  Rectangular  Panel  Card  (page  8,  Ref.  2) 

RH0CF  (E8.U,  Col.  55-62)  -  RH0CF  is  the  product  of  Yf  •>   the 

(between  RH0  and  Table     fluid  specific  weight,  and  c,  the 
^  '  sonic  velocity  of  the  fluid.  The 

c  -2      -1 

units  of  Yf>  are  lb„  -in.   -  sec. 

If  the  panel  is  not  in  contact 
with  the  fluid,  RH0CF  =  0. 
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IV.   SAMPLE  PROBLEM 

A  simply  supported  square  plate  is  subjected  to  a  step  pressure  load 

of  the  form 

p  =  o  t  =  o 

p  =  P  sin  2*  Sin  iM.  t  >  o  (22) 

a      a  —  x   ' 

Due  to  symmetry,  only  one  quarter  of  the  plate  is  considered:  The 

parameters  of  the  problem  are: 

E  =  10.U  x  10  psi  -  Young's  modulus 

Y  =  O.O965  Ib/inT  -  specific  weight  of  the  plate 

y  =  l/3  -  Pois son's  ratio 

h  =  0.1  in.  -  thickness 

a  =  20  in.  -  length  and  width 

P  =  0.01  lb/in? 

The  load  is  sufficiently  small  such  that  the  nonlinear  effects  are 

negligible.  The  plate  is  modeled  with  four  elements  as  shown  in  Fig.  1. 
The  equation  governing  the  damped  motion  of  the  plate  corresponding 

to  Eq.  16  is 

Dy^w  +  vH  w  +  pcw  =  Psin  tiL  Sin  iSL  (23) 

g  a  a 

where 

3  ,  k  h  k 

D  "        ,        27'      V  IT       2  ~~ 2~~ 2  k 

12(l-v    )  &x*  &&  3yH 

and  g  is  the  local  acceleration  due  to  gravity.     The  solution  to  Eq. 

23  is 
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w  =  wgt  [1  -  e"Ctut  cos  Cyi-Q     (l)t  +  <p)  /cos  cp}  (24) 

where 

tan  ^  =  -  C/  VI  -  C2 


*P  sin  H*  sin  ^ 

st   4DnU 


^   2Yh 


«  "  2   V  hV 
a 


V  hv 


when  the  plate  is  initially  at  rest.  The  displacement  at  the  center  of 
the  plate  given  by  Eq.  2k   is  plotted  in  Figures  2,  3?  and  k   as  a  function 
of  time  for  Q  =   0,  0.666  and  270  corresponding  to  zero  damping,  less  than 
critical  damping  and  very  heavy  damping  respectively.  The  corresponding 

values  of  gpc  for  the  fluid  are  0,  k,   and  1620  lb  /(in  -sec).  Also 
plotted  in  Figs.  2-k   are  the  results  from  BR-1HR.  The  input  data  sheets 
and  the  print  of  the  input  data  are  given  in  Fig.  5-   The  execution  time 
on  the  IBM  360/67,  FORTRAN  IV  -  Level  H,  was  8  min.  and  26  sec.  for  200 
time  steps  with  r  =  270.  The  run  with  damping  not  considered  took 
essentially  the  same  length  of  time. 
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V.   SUMMARY  AND  CONCLUSIONS 

The  finite  element  digital  computer  code  BR-1  developed  by  the  Northrop 
Corporation  for  predicting  the  effects  of  internal  air  blast  on  typical 
combat  aircraft  skin- rib- stringer  structures  has  been  modified  to  include 
the  effect  of  compressible  fluid- structure  interaction.  The  fluid- structure 
interaction  is  approximated  by  the  piston  theory  wherein  the  effect  of  the 
fluid  upon  the  structure  is  accounted  for  by  introducing  damping  to  the 
equations  of  motion  of  the  structure.  The  modified  code  is  called  BR-1HR. 
This  code,  in  conjunction  with  the  NWC  code  for  predicting  hydraulic  ram 
pressures,  can  be  used  to  predict  the  structural  response  of  aircraft  fuel 
tanks  subjected  to  penetrating  bullets  and  fragments. 

All  of  the  features  of  BR-1  exist  in  BR-1HR,  and  only  two  additional 
numbers  are  required  for  the  input  data.  The  modified  code  is  operational 
on  the  IBM  360/67  in  FORTRAN  IV,  level  H,  and  requires  290K  bytes  of 
storage.  A  sample  problem  was  executed  to  demonstrate  the  validity  of  the 
modified  code  for  zero  damping,  less  than  critical  damping,  and  very  heavy 
damping. 
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APPENDIX  -  A  STUDY  OF  THE  ACCUEACY  AND  STABILITY  OF  SEVERAL  NUMERICAL 
INTEGRATION  SCHEMES  FOR  THE  TRANSIENT  RESPONSE  OF  HEAVILY  DAMPED 
STRUCTURES 

Many  studies  have  "been  made  of  the  accuracy  and  stability  of 
numerical  integration  schemes  for  the  equations  of  motion  of  structural 
systems.  However,  most  of  these  studies  concentrate  on  the  response  of 
undamped,  or  lightly  damped,  systems.   Of  interest  here  is  the  response 
of  both  lightly  damped  and  heavily  damped  systems  since  both  kinds  of 
damping  can  occur  when  a  structure  is  vibrating  in  contact  with  a 
compressible  fluid  (6). 

The  equations  of  motion  of  the  system  under  consideration  are  given 
in  matrix  form  by  Eq.  16.  Three  different  finite  difference  schemes  for 
the  numerical  solution  of  these  equations  are  considered  here.  Only 
explicit,  non-iterative  schemes  are  considered  due  to  the  fact  that  the 
BR-1  code  uses  an  explicit  solution  procedure.  Two  of  the  three  schemes 
are  based  upon  a  two  variable  approach  using  {q*}  and  {v*},  where 

{4*}  =  {v*}  (B-l) 

Thus,  Eq.  16  can  be  given  in  the  form 

{>}  =  [M]_1  ({C}  -  [D]  {v*})  (B-2) 

First  Scheme 

Substituting  the  first  order  approximations  for  {v*}  and  {q*} 

{v*}t  =  ({v*}t    -  [v*}t  )  /(At)  (B-3a) 

i        i+1       i 
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{4*}t  =   ({q*}t  -   [q*}t   )/(At)   =  {Aq*}t     /(At)  (B-3b) 

i+1  i+1  i  i 

into  Eqs.   B-l  and  B-2  leads  to 

{v*}t         =  [I   -  M_1D   (At)]    {v»l        +   (At)    [M]"1  {CV  (B-I+a) 

i+1  i  i 

CA<l1+  =  At   {v*!  (B-l+b) 

Vl  Vl 

Eliminating  [v*]  from  Eqs.  B-k   results  in 

{Aq*}t    =  [I  "  M_1D  (At)]  {Aq*}.   +  (At)2  [M] _1  [C]        (B-5) 
i+1  i 

This  is  identical  to  the  scheme  used  in  BR-1  when  damping  is  not  considered, 
It  is  also  equivalent  to  the  scheme  where  the  acceleration  {q*}  is  approxi- 
mated by  the  conventional  central  difference  approximation,  Eq.  17a.  The 

two  variable  approach  given  by  Eqs.  B-k   may  be  more  desirable  than  Eq.  B-5 

2 
due  to  roundoff  error  considerations,  i.e.  (At)   in  Eq.  B-5  is  a  very 

small  number. 


Second  Scheme 

The  second  scheme  uses  Eq.  B-Ua  and  the  simple  forward  Euler  approxi- 
mation for  [v*]  in  Eq.  B-l 

{Aq*l+    =  (At)  {v*}  (B-6) 

i+l  % 

in  place  of  the  backward  approximation  of  Eq.  B-Ub.  The  solution  proce- 
dure is  to  compute  fv*},    using  Eq.  B-Ua  and  fAq*"L   using  Eq.  B-6. 

*i+i  Hi+i 

Third  Scheme 

The  third  scheme  uses  the  conventional  central  difference  approxi- 
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mations  for  both  [q*]  and  {q*} ,  i.e.  Eqs  .  15  and  17a.  This  gives  Eq.  18, 
reproduced  here  for  convenience 


{A<l*)t    =  [M  +  D  (At/2)]"1  ([M  -  D  fot/2)]  [Aq*lt  +  (At)2  {CI   )   (B-7) 
i+1  i  i 

This  scheme  is  also  identical  to  the  BR-1  scheme  when  damping  is  not 
considered.  A  two  variable  version  of  this  scheme  is 

{AV*}.    =  [M  +  D  (At/2)]  ([M  -  D  (At/2)]  {av*}   +  (At)  {C}   )     (B-8a) 


i+1 


t. 

i 


and 


[A<l*}t    =  (At)  {av*). 


\+l 


(B-8b) 


'i+1 


This  may  have  smaller  roundoff  error  than  Eq.  B-7  since  (^t)  has  been 
eliminated. 

The  Single  Degree  of  Freedom,  Damped  Oscillator 

The  equation  for  the  free  vibrations  of  the  single  degree  of  freedom, 
damped  oscillator  is 


mq  +  dq  +  hq  =  0 


(B-9) 


Applying  the  three  schemes  described  above  to  Eq.  B-9  leads  to 

qfc  -    (2   -  2   CtB  "  /)   <lt     +   (1  "  2Clo)qt  =  0 

i+1  i  i-1 


(B-lOa) 


(At)v 


'i+1 


1   "  2C 


to  "(0 


(At)v 


i 


(B-lOb) 


(1  +  Cm\         -  (2  -  »   )qt     +  (l  -  ^\         =  0  (B-10c) 

i+1  i  i-1 
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where 

=   (&t)u,       m  =  Vkfm        q  =  d/(2m(0) 


0) 


according  to  Eqs.  B-5,  B-Ua  and  B-6,  and  B-7,  respectively. 
The  solution  to  Eq.  B-9  can  he  given  in  the  form 


^     =  A  (e~&>  e^VS     ~  1\       =  (a^)1  (B-ll) 

where  A  is  an  arbitrary  constant  and  the  superscript  i  denotes  the  ith 
power.  The  solution  to  the  three  difference  schemes  for  an  arbitrary 
set  of  initial  conditions  can  be  obtained  by  assuming 

C4.  =  A\X  (B-12a) 

i 

(At)vt  =  B\X  (B-12b) 

i 

where  \,  A  and  B  are  unknown  constants.  Substituting  Eqs.  B-12  into 
Eqs.  B-10  and  solving  for  \   for  each  scheme  lead  to 

\±  =   1  "  C5  "  To2/2  t  5  \/(C  +  (B/2)2  -l  (B-13a) 


X2  =  1  "  C(S  T  5  \/c2  "  1  (B-13b) 


X3  =  (1  -  u52  /2  -  5  \/"c2  +  ^A  -  1  )  /(I  +  C£)        (B-13c) 

When  the  discriminant  in  Eqs.  B-ll  and  B-13  is  positive  the  solution 
consists  of  damped  motion  only.  When  it  is  negative  the  motion  is 
damped  and  oscillatory.  Thus,  a  zero  discriminant  defines  the  limit 
of  the  oscillatory  behavior. 
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The  accuracy  of  the  three  numerical  schemes  can  be  evaluated  by- 
comparing  \    ,  \  and  \     with  \     for  several  values  of  g  and  ^.   The 

values  of  the  ratios  of  the  numerical  solution  to  \     are  given  in  Table 

c 

B-l  for  £  -   0,  0.5,  5,  and  500  and  ^  =  0.1.  This  value  of  (j  corresponds 

to  a  time  step  of  1  ^  sec  when  ^  =  100,000  rad/sec  or  a  time  step  of  1 

msec  when  (0  =  100  rad/sec,  etc.,  i.e.  the  solution  is  computed  ten  times 

over  the  time  interval  l/(0  or  20^  times  over  the  undamped  natural  period 

2tt/u>*  ^e  closer  the  ratio  in  Table  B-l  is  to  one,  the  closer  the 

numerical  eigenvalue  is  to  the  correct  eigenvalue. 

The  numerical  stability  of  each  scheme  can  also  be  determined  from 

Eqs.  B-13.  When  |\|>1  the  numerical  solution  will  be  unstable.  The 

upper  limit  on  £  for  stability  can  be  determined  for  a  given  value  of  Q 

by  equating  |\|  to  one  and  solving  for  ty.  When  the  discriminant  is 

positive  \  =  -1  is  the  limiting  value  and  the  negative  sign  in  the  =f 

/  2    2 

applies.  When  the  discriminant  is  negative  |\|  =Jx     +   y  where  x  and 

y  are  the  real  and  imaginary  parts  respectfully.  The  results  for  the 
limiting  values  of  £  for  an  oscillatory  solution  and  for  numerical 
stability  are  given  in  Table  B-2. 
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Scheme 

Oscillatory  limit 

Stability  limit 

*1 

£  =  2(1-C) 

d5  =  2(\/c2  +l'-C> 

x2 

c  =  1 

m  =  2£,  c  *  x 

n  =  2c-  2  v  c2-1  » 

x3 

£  =  2 

5  =  2  \/l  -c2 

TABLE  B-2  Limits  on  £  for  an  Oscillatory  Solution  and  a  Stable 

Solution 


Note  that  \  is  unstable  for  any  non-zero  value  of  $  when  £=0  and  that 
2  is  the  maximum  limit  on  fc   for  all  three  schemes.  Also  note  that  the 
limit  on  \_  is  the  same  as  that  on  the  BR-1  routine,  even  with  damping, 
and  that  this  is  the  least  restrictive  scheme.   Consequently,  based 
upon  these  accuracy  and  stability  considerations  the  third  scheme  is 
selected. 
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